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We present a mathematical characterization of hyperbolic gauge pathologies in electrodynamics 
and general relativity. We show analytically how non-linear gauge terms can produce a blow-up of 
some fields along characteristics. We expect similar phenomena to appear in any other gauge field 
theory. We also present numerical simulations where such blow-ups develop and show how they 
can be properly identified by performing a convergence analysis. We stress the importance of these 
results for the particular case of numerical relativity, where we offer some cures based on the use of 
non-hyperbolic gauges. 
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During the last decades, gauge field theories have be- 
come the paradigm of fundamental physics. In such the- 
ories, gauge invariance implies that the 'physics' is inde- 
pendent of our choice of gauge. Yet in some cases, of 
which general relativity (GR) is a very good example, it 
can be difficult to separate the gauge from the physical 
degrees of freedom. This is of particular importance in 
numerical simulations of non-linear field theories where 
one needs to distinguish physical singularities from pure 
gauge pathologies. In such cases, exact solutions are not 
available, and the simple idea of a change of gauge can 
become a daunting task. Because of this, understand- 
ing the properties of particular gauge choices becomes 
an issue of great relevance. 

In a recent paper it was shown, to the surprise of 
many in the numerical relativity community, that cer- 
tain gauge choices (throughout we use the term 'gauge' 
applied to relativity to mean only slicing conditions) in 
relatively simple scenarios could lead to the development 
of what were there called 'coordinate shocks'. Here we 
will show that the term 'shock' used in that paper is mis- 
leading and will refer to them instead as gauge patholo- 
gies. The recent Bona-Masso hyperbolic formulation of 
the Einstein Equations (|] allowed the author of reference 
to study the structure of the gauge condition and how 
its non-linearity could cause gauge pathologies. It was 
shown numerically that these pathologies seemed to oc- 
cur, but no proof that they were real discontinuities was 
given. This work left some basic questions unanswered: 
How can we characterize these pathologies mathemati- 
cally and numerically? How generic arc they? Do similar 
phenomena occur in other gauge theories? 

Here we give for the first time a mathematical char- 
acterization of hyperbolic gauge pathologies (where by 
hyperbolic gauge we mean choices where the gauge is 
evolved using a hyperbolic equation) based on the the- 
ory of nonlinear waves H . We focus our attention on two 



cases: electrodynamics (ED) and spherically symmetric 
GR. We analyze the structure of the equations and show 
how non-linear gauge terms can produce a blow-up of 
some fields along characteristics ||. Such a blow-up in- 
dicates a gauge singularity and corresponds to the 'coor- 
dinate shocks' of reference §. The blow -ups described 
here are much stronger singularities than shocks since 
some of the propagating fields become in fact infinite. 
We also show how the blow-ups can be identified numer- 
ically by analyzing the convergence properties of compu- 
tational simulations. Finally, we propose some "cures" 
for these pathologies based on changing the hyperbolic 
nature of the gauge condition. Thought we only study 
the cases of ED and GR, we expect similar phenomena 
to appear in any other gauge field theory. 

Electrodynamics — Surprisingly, gauge effects can pro- 
duce blow-ups along characteristics even in simple sys- 
tems such as ED. We are not aware of any previous anal- 
ysis of gauge pathologies in ED, probably because well 
behaved gauge choices are intuitively clear there. In the 
following, we will follow closely the notation of Refs. jjJU ■ 

We will write the equations of ED as a first order initial 
value problem in the following way: 

d t A i = -(E i +^ i ) , (la) 
d t D ij = -d i {E j +i) j ) , (lb) 
d t E i = d i trD-J2 j 9 j D ji , (lc) 

where cf> and A are the scalar and vector potentials, E the 
electric field, and where we have introduced the quan- 
tities Dij := diAj and ipi := di<f>. We also have one 
constraint which in vacuum takes the form V • E = 0. 

We are now free to decide how {4>, ipi} will evolve, i.e. 
we are free to choose the gauge. Here we make a choice 
similar to the one we will make later in relativity 

d t cf> =-/(<£) txD, (2) 
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with /(</>) > but otherwise arbitrary (/ = 1 corre- 
sponds to the familiar Lorentz gauge). It is through the 
gauge function / that we introduce a non-linearity into 
ED, which is otherwise a completely linear theory. 

The resulting system of equations can be shown to be 
hyperbolic (in the sense of having a complete family of 
eigenfields). For a given direction x we find the charac- 
teristic structure: 4 fields propagate along the physical 
light cones (speed ±1) 

w p± := (E p + ip p ) ± D xp , p jt x , (3) 

and 2 fields propagate with the 'gauge speeds' ± \fj 

uj 9± :=^ x ±^/ftxD . (4) 

The remaining 13 fields move along the time lines. 
Let us now define Q g± '■= \[j We find 

d t n g± ±^]d x n g± = ±-£j(n g l + n g+ £i g J s ) (5) 

For /' 7^ 0, the quadratic term in S! s± can produce a 
blow-up along a characteristic in a finite time. We can 
predict the time when this will happen if we restrict our- 
selves to the case £l g _ = 0, and take / = e a ^ . Along a 
characteristic we will have 



9+ ' 



dt 4 

which can be easily integrated to find 



9 + 1 - ail t/4 



(0) 



(7) 



For a^o > 0, a gauge pathology (a blow-up of Sl g .) 
will appear at a finite time T given by 



T = A/aQo 



(8) 



We have constructed a numerical code to evolve this 
system of equations and reproduce this blow-up. Fig. |l| 
shows the results of a one-dimensional simulation for 
which a shock is expected at T = 19.92. The first two 
plots show the initial (dotted line) and final (solid line) 
values at t — 20 of <f> an d ip x - A sharp gradient has 
developed in 4>, while a large spike has appeared in ip x . 

Now, if we had not performed the mathematical anal- 
ysis of the system, we could ask ourselves how can we 
know that these features correspond to a real blow-up 
and not to the development of large but smooth gradi- 
ents. After all, nonlinear systems may develop features 
that are difficult to resolve numerically. As an exam- 
ple, one of the key successes of numerical relativity has 
been the discovery of critical phenomena j| . Confidence 
in these results required advanced adaptive mesh refine- 
ment and the careful use of convergence tests ||. Here 
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FIG. 1. Numerical simulation of a gauge wave in electro- 
dynamics. The first two plots show the initial (dotted line) 
and final (solid line) values at t — 20 of the scalar potential (f> 
and its derivative ip x . The lower plots show the convergence 
rates as a function of time. 



we have performed similar convergence studies: we con- 
sider a series of different resolutions and find the rate at 
which a global measure of their relative errors converges 
to zero (the 'self-consistent' convergence rate). Our code 
uses a second order accurate scheme, so we expect the 
convergence rate of our solutions to be 2. The lower two 
plots of Fig.|l| show the convergence rates for </> and ip x 
obtained from runs at three resolutions (5000, 10000 and 
20000 points). As expected, for t < 20 the convergence 
rate is close to 2, while after that it drops dramatically 
indicating loss of convergence. Moreover, the spikes in 
ip x become larger and larger with increased resolution. 
We then conclude that they correspond to a real infin- 
ity and not just a large gradient. The lesson is clear: we 
can characterize numerically the appearance of a blow-up 
using convergence analysis. 

It is important to consider the effect that the choice 
of finite difference method has on the results described 
above. When we used standard methods like a leap frog 
scheme, the simulations crashed very soon after the spike 
in i/j x started to develop. Using instead shock capturing 
Q techniques we were able to follow the evolution much 
further, which allowed us to determine more precisely the 
time of the blow-up. Notice that this further evolution 
is non-physical: the shock capturing techniques help to 
maintain the simulation stable but this is of no physical 
relevance since we the true solution is not really a shock 
wave. A real infinity has developed and no numerical 
method will converge after that. 

Spherically symmetric general relativity- Let us now 
consider the case of GR in spherical symmetry. The ba- 
sic variables are the lapse function a, the spatial metric 
components {g rr ,gee} and the extrinsic curvature com- 
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ponents {K rr , Kgg}. Since we are interested in hyper- 
bolic gauges we will use the following gauge condition 0] 



d t a = 



] f(a)trK 



(9) 



with f(a) > but otherwise arbitrary (/ = 1 now corre- 
sponds to harmonic slicing). It can be shown that with 
this gauge condition the evolution equations of GR can 
be written in first order hyperbolic form 0] . 

The particular form of the equations for spherical sym- 
metry and its characteristic structure can be found in [Q . 
There it was shown that there are two families of travel- 
ling modes, one that moves with the speed of light, and 
one that moves with a 'gauge speed' that depends on the 
value of the function / and that reduces to the speed 
of light for harmonic slicing. In general, both types of 
travelling modes can develop gauge pathologies. Here we 
will concentrate in the particular case of harmonic slicing 
(/ = 1) where only one type of pathology appears. Note 
that this gauge is the preferred choice for many of the 
new hyperbolic formulations of the Einstein equations 
The analysis of the system for other forms of / will 



be considered elsewhere 12 



We will start by defining 



fi± := a /gee i^Kgg ± D r ee/g. 



1/2 
rr 



(10) 



where D r ee '.= 1/2 d r ggg. ^From the form of evolution 
equations given in Q it is easy to show that, for / = 1, 
we will have 



d t n± ± a/gli 2 d r n± = - fi± + "Vffe 



(11) 



We see that the 17 ± represent outgoing (+) and ingoing 
(-) modes travelling with the speed of light. Unfortu- 
nately, the last term on the right hand side of the above 
equations makes it impossible to separate them, and so 
prevents one from predicting when a blow-up caused by 
the quadratic source terms might occur. Nevertheless, 
one can see that the quadratic term in 17± appearing 
in ([ll]) is only dangerous for negative values of fl± . More- 
over, it is easy to see that if 17± is initially positive it 
will not change sign. Now, if our initial data is time- 
symmetric (Kij — 0) and such that the metric function 
ggg is monotonic, then we will have ± Cl± > 0. This 
means that no blow-ups can develop for outgoing modes: 
Only ingoing modes can produce a gauge pathology. Of 
course, whether they will or not depends on the precise 
form of the initial data. 

In reference [[t| it was shown that pathologies appeared 
for a particular choice of initial data in a black hole space- 
time. Here we give another example of how pathologies 
can form from apparently simple initial data. We choose 
the standard initial data for the metric and extrinsic cur- 
vature that has been successfully used for most black hole 
simulations to date P-[To| . It corresponds to an isotrop- 
ically sliced Schwarzschild black hole with time symme- 
try. The key difference is that, instead of choosing an 



initial lapse that satisfies the maximal slicing condition, 
we choose the following 'Gaussian' profile: 



a = 1 - A exp [((r - r )/aY 



(12) 



We have used 3 independent numerical codes to evolve 
the system, one based on the standard Arnowitt-Deser- 
Misner (ADM) formulation jll|] that uses a simple leap 
frog scheme, and two based on the Bona-Masso hyper- 
bolic formulation || that use shock capturing methods 
[^J. All three codes produce similar results with one im- 
portant difference: the ADM code crashes soon after the 
pathologies start to develop, while the other two codes 
are capable of continuing past this point. 

We report the results of a particular simulation ob- 
tained using one of the hyperbolic codes. The first plot 
in Fig.|] shows the initial (dotted line) and final (solid 
line) values of the lapse at time t = 15M for the case 
{A = l,ro = M/2,<7 = 6M,p = 2}, with M the mass 
of the black hole. The next plot shows the same for the 
conformal metric function g rr . Notice how both the lapse 
and the radial metric develop large spikes, the sizes of 
which increase with resolution. The bottom plots show 
the global convergence rate of the lapse and the hamil- 
tonian constraint obtained from runs at 4000, 8000 and 
16000 grid points, with the outer boundary located at 
40M. Clearly, we have a gauge pathology, whose blow- 
up time seems to be t = (14 ± 1)M. 
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FIG. 2. Numerical simulation of a spherically symmetric 
black hole spacetime. The first two plots show the initial 
(dotted line) and final (solid line) values at t = 15M of the 
lapse a and the conformal metric function g rr . The lower 
plots show the convergence rates for the lapse and the hamil- 
tonian constraint. 

Although the simulation presented here corresponds 
to harmonic slicing, gauge pathologies also develop with 
other hyperbolic gauge choices. An important difference 
is that in the harmonic case the lapse becomes infinite, 
which indicates that the time slicing becomes null. In the 
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other cases, the lapse becomes discontinuous but remains 
finite, and the time slicing develops a kink instead. We 
have performed many simulations studying the parame- 
ter space {A, a, p} and found similar results. Details of 
all these studies will be given elsewhere Q. 

Another crucial aspect of the problem is that of the 
numerical resolution of the gauge pathologies. We have 
found that if we evolve the previous system with only 
200 grid points, the pathologies do not seem to form. 
The lapse grows until a certain value is reached and then 
it propagates out in a smeared manner due to the large 
numerical viscosity. The solution "looks good" , although 
it is non-physical, as a proper convergence test reveals. 

The impact of these results on three dimensional (3D) 
numerical relativity should not be underestimated. Note 
that in the previous examples we have used thousands 
of points to be able to show very sharp fall-offs of the 
convergence rate and have a good estimate of the blow- 
up time. Even if this sharpness will not be possible at 
the resolutions currently available for 3D computations, 
at medium and low resolutions one can already see that 
the convergence fails at late times. 

We should stress again the fact that the development 
of these pathologies depends crucially on the form of the 
initial data. For different choices of the initial lapse func- 
tion, one can find that harmonic slicing is perfectly well 
behaved, as the simulations presented in reference @] 
show. In fact, one can even find explicitcly a harmonic 
slicing of a black hole in which the metric is static jl3) . 

"Cures"- The appearance of gauge pathologies might 
seem to put into question the practical value of non-linear 
hyperbolic gauges in numerical studies of gauge field the- 
ories. After all, in a general situation it might be very 
difficult to know a priori if our initial data will develop 
such a pathology. 

For ED, the solution is clear: use a gauge that decou- 
ples the characteristic speeds from the dynamics, i.e. use 
the Lorentz gauge. Unfortunately, this will not work in 
relativity where the characteristic speeds cannot be de- 
coupled from the dynamics. We can think of at least two 
different ways to solve the problem. Both involve chang- 
ing the character of the equations for the gauge only: all 
other equations (the 'physics') remain hyperbolic. 

The first approach implies using an elliptic gauge con- 
dition of which maximal slicing is the best known exam- 
ple. One can then either use an elliptic gauge always or, 
in cases where it might be of interest, use a hyperbolic 
gauge for some time and then switch to an elliptic gauge 
when a pathology is about to form. We have tested this 
idea and found that it works very well in practice ]l2] | . 

The second approach consists of adding dissipation to 
our gauge condition. We will then have a parabolic equa- 
tion and intuition tells us that this should prevent the 
pathologies. We then propose the gauge condition 

d t a = -a 2 [f(a)trK- £(a) V 2 a] , (13) 



with /, £ > but otherwise arbitrary. Notice that the 
numerical treatment of the diffusion term requires either 
a very stringent Courant condition or the use of implicit 
techniques. Note also that this term should be kept with 
the same coefficient for all resolutions of a given simu- 
lation, as it does not correspond to a simple "artificial 
viscosity" ]l0[ ], but rather to an explicit change of the 
character of the gauge condition. We have tried this con- 
dition in spherically symmetric GR and found that it also 
prevents the development of gauge pathologies |Q. 

In conclusion, we have presented for the first time a 
characterization of hyperbolic gauge pathologies in ED 
and GR. We have shown how the coupling of character- 
istic gauge speeds to the dynamics produces a nonlin- 
ear blow-up mechanism, and how a careful convergence 
analysis can indicate the appearance of such a blow-up 
in numerical simulations. The origin of these pathologies 
is in the finite speed of propagation of the gauge modes, 
and therefore a way to avoid them is the use of elliptic 
or parabolic gauges with infinite gauge speeds. 

Thought we have concentrated in the cases of ED and 
GR, we expect similar pathologies to arise in any other 
gauge field theory. Because of this we feel that further 
mathematical study of these phenomena will be of fun- 
damental importance for future numerical simulations of 
non-linear field theories. 
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